/* THIS PROGRAM IDENTIFIES THE MEMBERS OF A FAMILY WHICH ARE ASSOCIATED TO
   A BODY *idesig_to_start* AT A LEVEL *limit* (TO BE SET BELOW) 
   WORKS ON ANALYTICAL ELEMENTS OF MILANI AND KNEZEVIC, HCM OF ZAPPALA ET AL. */

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#include "nrutil.h"
#include "nrutil.c"

#define PI2 6.28318530717959 
#define AU 1.49597870e11 // IN METERS

#define GAUSS 0.01720209895
#define GM (GAUSS*GAUSS*365.25*365.25) //UNITS: JULIAN YEARS, AU, SUN MASS

#define FACTOR (AU/365.25/86400.) // CONVERSION FACTOR FROM AU/YR TO M/S

#define MAXN_IN_LIMIT 1000 // MAXIMUM NUMBER OF BODIES ALLOWED WITHIN A *LIMIT* ABOUT ANY BODY 
                           // IF EXCEEDED, PROGRAM GIVES A RUN-TIME ERROR MESSAGE AND EXISTS  

static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

void fgetfield(FILE *fp,char s[],int length);

int main(void){

  char filename[1000];
  char field[500];
  int RFL,QCM,QCO;
  long int i,j,k,l,inew,ide,ntpnumb,ntpmult,ntp;
  long int *idesig,*nassoc,**assoc,nfamily,*family;
  long int idesig_to_start, k_to_start;
  double *sema,*ecc,*sinc,g,s,*mag;
  double limit,sqrfac,metrics,maglimit;
  FILE *fp;

  //********************************************************************************

  printf("Enter input filename:  ");
  scanf("%s",filename);
  printf("Enter id of the central body (e.g. 490 is Veritas) : ");
  scanf("%ld",&idesig_to_start);
    
  //  idesig_to_start = 158; // ENTER HERE THE DESIGNATION NUMBER OF THE BODY YOU WISH 
                         // TO ASSOCIATE THE FAMILY WITH 

  printf("Enter cut-off : ");
  scanf("%lg",&limit);

  printf("Enter limiting magnitude of bodies to be considered: ");
  scanf("%lg",&maglimit);

  // limit = 100; // SEPARATION LIMIT, IN M/S, MUST BE LARGER THAN QRL (ZAPPALA ET AL.)
  
  limit = limit*limit;

  //********************************************************************************

  sqrfac = SQR(FACTOR); // TO AVOID SQUARE ROOT IN METRICS, WORK WITH METRICS^2  

  ntpnumb = 23845;
  ntpmult = 42244;

  ntp = ntpnumb + ntpmult;

  ntp=0;

  if ((fp = fopen(filename, "r"))== NULL ) 
    {
      printf("E 00 File not found\n");
      return 0; //if does not succed, return 0
    }
  
  // count number of lines in the file: only ephemerides must be present
  while (1)
    {
      if (NULL==fgets(field,256,fp)) break;
      if (field[0]=='#') continue;
      ntp++;
    }

  fclose(fp);//close the file in order to rewind the pointer
  
  sema = dvector(1,ntp); ecc = dvector(1,ntp); sinc = dvector(1,ntp);
  mag =  dvector(1,ntp);
  idesig = lvector(1,ntp);

  nassoc = lvector(1,ntp);
  assoc = lmatrix(1,ntp,1,MAXN_IN_LIMIT);
  family = lvector(1,ntp);

  // READ THE PROPER ELEMENTS OF NUMBERED ASTEROIDS 
  if ((fp = fopen(filename, "r"))== NULL ) 
    {
      printf("E 00 File not found\n");
      return 0; //if does not succed, return 0
    }

  for(i=1;i<=ntp;i++){
       if (NULL==fgets(field,256,fp)) break;

    sscanf(field,"%ld %lf %lf %lf",
    	   &idesig[i], &sema[i],&ecc[i],&sinc[i]);
    mag[i]=20;
  }

  fclose(fp);

  // WHAT IS THE BODY TO START WITH ? 
  i = 1;
  while(idesig[i++] != idesig_to_start && i<=ntp);

  if(i>ntp){ 
    fprintf(stderr,"The central body, %ld, not found ...", idesig_to_start);
    fprintf(stderr,"Exiting now ...\n");
    return 0;
  }

  k_to_start = i-1;

  // NASSOC[I] IS # OF BODIES WITHIN *LIMIT* DISTANCE TO BODY *I*, INITIALLY SET TO ZERO
  for(i=1;i<=ntp;i++)
    nassoc[i] = 0; 
  
  // SELECTION OF BODIES LINKED WITH *idesig_to_start*
  k = k_to_start;         
  nfamily = 1;
  family[nfamily] = k;       // add *idesig_to_start* as 1st body to the family list 

  i = 1;  /* run with *i* over the current list of family members, 
	     attaching additional branches */
  
  while(i <= nfamily){  // stop if *i* is larger than the current family list: no more branches
   
    k = family[i];   /* run over the bodies associated to *i* body */
          
    for(j=1;j<=ntp;j++)
      
      if(j != k){
	  
	// D1 METRICS SQUARED (!) OF ZAPPALA ET AL. 1994
	metrics = 
	  1.25*SQR( 2.00*(sema[j] - sema[k])/(sema[k] + sema[j]) ) + 
	  2.00*SQR(  ecc[j]  - ecc[k]           ) +
	  2.00*SQR(  sinc[j] - sinc[k]          );  	  
	  
	  metrics = 2. * metrics*GM/(sema[j] + sema[k]) * sqrfac;
	  
	  // IF J-TH BODY IS CLOSER THAN *LIMIT*, ASSOCIATE IT TO I-TH BODY
	  if(metrics <= limit && mag[j] <= maglimit){
	    nassoc[k] ++;
	    
	    if(nassoc[k] > MAXN_IN_LIMIT){
	      fprintf(stderr,"Number of associated bodies to %ld exceeds MAXN_IN_LIMIT ...\n",k);
	      fprintf(stderr,"Increase MAXN_IN_LIMIT, Exiting ...");
	      return 0;
	    }
	    
	    assoc[k][nassoc[k]] = j;
	  } 
	
      }      
        
    for(j=1;j<=nassoc[k];j++){   
      
      inew = 1;      /* see if bodies associated to *k* are already in the family list */
      for(l=1;l<=nfamily;l++)
	if(assoc[k][j] == family[l] )
	  inew = 0;

      if(inew == 1){    /* if this is a new body, add it to the family list */
	nfamily++;
	//	fprintf(stderr,"%ld\n",nfamily);
	family[nfamily] = assoc[k][j];
      }
  
    }

    i++;
  }

  // OUTPUT THE FAMILY LIST 
  fp = fopen("family.list","w");
  for(i=1;i<=nfamily;i++){

    ide = family[i];

    fprintf(fp,"%ld %ld %lg %lg %lg\n",
	   i,idesig[ide],sema[ide],ecc[ide],sinc[ide]); 

  }
  fclose(fp);

  return 1;
}

void fgetfield(FILE *fp,char s[],int length){
  int i;
  
  for(i=0;i<length;i++)    
    s[i]=fgetc(fp);
  s[i]='\0';
}    
